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This  report  considers  spectral  estimation  and  extrapolation  techniques 
for  discrete  time,  band  limited  signals,  (l.e.,  signals  whose  bandwidth  Is 
less  than  ^ cycles;  if  T sec.  is  the  sampling  Interval)  which  are  observable 
only  for  a Tinite  duration.  The  objective  Is  to  determine  the  spectrum 
(or  power  spectrimi)  of  these  signals.  It  Is  shown  that  the  estimated  spectrum 
can  be  Improved  considerably  (over  a perlodogram  or  Maximum  entropy  spectriim) 
by  first  extrapolating  the  given  observations  beyond  the  observation  Interval. 
Also,  we  consider  the  problem  of  extrapolation  of  signal  In  the  presence  of 
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noise  or  other  interfering  signals. 

Several  new  results  anu  algorithms  are  presented.  First,  it  is  shown 
that  some  of  the  existing  extrapolation  methods  for  continuous  signals  when 
extended  to  sampled  data  do  not  converge  to  the  exact  original  time-unlimited 
signal.  Rather,  one  only  expects  to  get  a minimum  norm  least  squares  estimate. 
And,  we  find  that  Papoulis'[8]  iterative  extrapolation  algorithm  Is  a special 
case  of  a gradient  algorithm  with  linear  convergence.  It  Is  shown  that  an 
infinite  extrapolation  matrix  Introduced  in  [10]  does  not  exist  and  is  ill- 
conditioned  at  best  when  approximated  to  a finite  matrix.  The  new  extrapola- 
tion algorithms  Include  a discrete  prolate  spheroidal  wave  function  (PSWF) 
expansion,  a conjugate  gradient  Iterative  algorithm,  a mean  square  extrapola- 
tion filter  and  a recursive  Kalman  filter  type  extrapolator . . The  latter  two 
algorithms  also  consider  the  noise  statistics  In  extrapolatim  of  the  signal. 
Several  examples  are  given  and  comparisons  are  made.  \ 
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I.  INTRODUCTION 

Spectral  estimation  refers  to  the  problem  of  estimating  the  spectral 
density  function  of  a stationary  random  signal  which  is  observable  only  . 
over  a finite  duration.  For  a deterministic  signal  it  implies  estimation 
of  its  magnitude  spectrum.  In  either  case,  if  the  signal  were  known  over 
the  infinite  interval,  the  Fourier  transform  of  the  signal  or  its  auto- 
correlation would  immediately  yield  the  spectrum.  Thus,  any  estimated  spectrum 
is  equivalent  to  specifying  the  signal  or  its  autocorrelation  outside  the 
observation  interval--i .e. , its  extrapolation. 

In  this  report  we  consider  several  algorithms  for  extrapolation  and 
spectral  estimation  of  discrete  time  signals.  First,  we' briefly  review 
the  maximum  entropy  (ME)  or  the  linear  predictive  autoregressive  (AR) 
method,  and  some  iterative  and  matrix  inverse  based  extrapolation 
algorithms  developed  recently  by  Papoulis  [8],  Sabri  and  Steenaart  [10], 
and  Cadzow  [11]. 

The  new  results  presente^^  here  are  as  follows. 

1)  Papoulis'  iterative  algorithm  applied  in  discrete  time  domain  converges 
to  an  extrapolated  signal  which  is  a minimum  norm  least  squares 

type  solution.  It  is  seen  to  be  a special  case  of  a one  step  gradient 
algorithm,  and  has  linear  convergence.  The  convergence  of  this  can 
be  improved  by  suitably  modifying  it  to  a steepest  descent  algorithm. 

2)  Sabri  and  Steenaart  [10]  have  reformulated  Papoulis'  iterative  algorithm 

in  terms  of  an  extrapolation  matrix  operator  which  yields  the 
extrapolated  signal  when  it  operates  on  the  given  time  truncated 
signal.  It  is  proven  that  the  infinite  operator,  E^  does  not  exist, 
but  its  finite  truncation  E^  exists,  but  it  is  ill-conditioned. 
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3)  It  is  known  that  a continuous  (time)  band  limited  signal  given  over 

a finite  observation  interval  can  be  extrapolated  exactly  outside  this 
interval  by  means  of  Prolate  Spheroidal  Wave  Functions  (psWF).  We 
show  that  for  the  discrete  time  case  a similar  expansion  arises  when 
we  consider  the  minimum  norm  least  squares  extrapolation. 

Then  we  present  three  other  algorithms  which  are  as  follows: 

4)  Conjugate  Gradient  Iterative  Extrapolation 

5)  Minimum  Mean  Square  Extrapolation  Filter 

6)  Recursive  Least  Squares  Extrapolation  Filter 

The  conjugate  gradient  method  is  an  iterative  algorithm  which  yields 
a psuedo  inverse  extrapolation  operator.  Compared  to  the  earlier  iterative 
methods  [8-11],  this  algorithm  converges  quite  rapidly.  The  minimum  mean 
square  extrapolation  algorithm  is  designed  for  applications  where  the 
observed  band  limited  signal  is  contaminated  by  wideband  white  noise.  It 
yields  a simple,  Wiener  filter  type,  extrapolation  operator  which  requires 
inversion  of  a matrix  whose  size  is  equal  to  the  number  of  samples  in 
the  observed  signal.  No  iterations  are  required  here  and  the  algorithm 
is  shown  to  reduce  to  the  matrix  inverse  algorithm  of  Cadzow  [11] 
as  the  additive  noise  power  goes  to  zero.  Finally,  the  recursive 
least  squares  algorithm  is  a Kalman  filter  based  method  where  the 
extrapolated  signal  estimate  is  updated  recursively  as  a new  observation 
sample  arrives.  The  latter  two  methods  are  applicable  in  the  presence 
of  noise  and  yield  stable  results.  Finally  these  algorithms  are  shown  to 
be  applicable  to  problems  where  one  needs  to  discriminate  as  well  as 
extrapolate  an  interfering  signal  and  a desired  signal. 

Several  examples  are  considered  to  compare  the  various  algorithms. 


A 


Tirffii~ni~iitfi''i  • 


II.  THE  MAXIMUM  ENTROPY  METHOD  [i!-7] 

Let  {u^.}  denote  a real,  zero  mean,  stationary,  Gaussian  random 
process  whose  covariance  function  is  defined  as 


We  know  r only  on  a finite  window  W defined  as 
m 

W = {-p<m<j)}  . (2) 

The  maximum  entropy  method  extrapolates  r^^  outside  W by  maximizing  the 


entropy 


e ^ ilnS{f)df 

-1/2 


under  the  constraint 


S(f)ej2iTmfdf  ^ 


The  solution  gives  the  maximum  entropy  spectrum  as 


S(f)  ■ 


This  could  be  written  as 


^ 

2 -jZufm  ^ 
) a e 

m=0 


where  the  a^^  and  are  related  by 


3 n, 

-m  m 


min(p-m,m) 

I 

k=max[0,-m] 


“k+m“k  • 


The  coefficients  (a^)  are  determined  by  solving  (4)  and  (6)  which  is  equivalent 
to  solving  (9)  below.  Alternatively,  the  fu^}  could  be  characterized  as  an  AR  process 

‘^k  = s“nVn  " ^k 
n=l 


(8) 
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where  6^  = By  writing  the  Yule-Walker  equations  for  (8)  it  is  a 

simple  matter  to  show  that  are  obtained  by  solving  a set  of  simultaneous 
linear  equations* 

Ra  = -6^I  . 6^  = l/(R’bij  (9) 


where  a and  1^  are  (p+1)  x 1 vectors  and  R is  a (p+1)  x (p+1)  covariance 
matrix  with  entries  corresponding  to  covariances  on  the  window  W,  i.e.. 


*For  a positive  definite  matrix  R,  are  guaranteed  to  be  such  that  S(f) 
is  positive  and  (8)  is  asymptotically,  a stationary  random  process. 
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Strictly  speaking,  the  covariance  values  m e W)  should  be  known 
exactly.  In  practice,  one  only  knows  data  values  on  a finite  window. 

Then,  the  covariances  could  be  estimated  as* 

1 M-|m| 

’ -B- “kVlml-  '’O’ 

where  M is  the  size  of  the  data  window.  For  large  M » p,  reasonable 
estimates  of  {r^^ could  be  expected. 

Note  that  this  method  does  not  require  {u^. } to  be  bandlimited  (with 
respect  to  the  Nyquist  rate).  Also,  the  spectral  density  function  is,  in 
view  of  (5)  and  (6),  an  all  pole  model.  Thus,  if  the  given  observations 
were  of  the  form 

yk=Uk*nk  (11) 

where  is  a white  noise  process  or  another  signal  (e.g.  clutter  noise 
which  could  be  modeled  by  an  AR  process),  the  spectrum  of  {yj^}  would  not 

be  an  all  pole  model  and  may  have  to  be  approximated  by  a very  high 
order  all  pole  model . 

Example  1(a):  Although  there  are  many  examples  where  the  maximum  entropy 
method  could  be  applied  successfully  [3,5,7]  we  consider  a case  where 
it  does  not.  We  assume  the  observations  to  be  given  by 

yk  = Sk  + Ck  * n^  (12) 

where  Sj^  represents  a bandlimited  signal  whose  spectrum  lies  in  the 
interval  [f2.f3]  and  [-f2.-^3]  and  is  an  interference  signal  band- 
limited  in  the  interval  [-fi,fi]  and  nj^  is  a white  noise  process. 

*In  estimating  r^^^,  the  divisor  of  M,  rather  than  M-jm]  is  recormended . 
Although  this  results  in  a biased  estimate  of  r^^,  it  yields  a positive 
definite  sequence  {rj^}  so  that  R is  positive  definite  and  the  resulting 
spectra  is  positive,  b^e  Parzen  [15]  for  details. 
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Figure  1(a)  shows  the  spectra  of  the  various  signals.  Figure  1(b)  and  1(c) 
show  the  interference  signal  C|^  and  the  actual  signal  Sj^,  modeled  as 

S|^  = 1 .69  Sin( . 39TTk) . 

Figure  1(d)  shows  the  17  samples  of  the  available  observations. 

The  signal  to  interference  signal  (to  be  called  clutter)  ratio,  which 
is  defined  as 


cro*  1 Peak  to  Peak  Value  of  Signal 
bLK  - log^g  clutter  ’ 


(13) 


is  -4.1  dB  and  the  signal  to  noise  ratio,  SNR,  defined  similarly  is  19  dB. 
Figure  1(e)  shows  the  maximum  entropy  spectrum  estimate.  A peak  is 
expected  at  the  position  marked  by  the  arrow.  At  this  point  the  signal 

estimate  is  30  dB  below  the  clutter  peak  and  is  indistinguishable  from 

o 

the  interfering  signal.  Figure  1(f)  shows  the  spectrum  estimated  by 
directly  evaluating  the  Fourier  spectrum  (i.e.  the  periodogram)  as 


8 

I y.exp(-j2Trfk) 
k=-8 


2 


-i<  f <i 
2 - - 2 


(14) 


Equation  (14)  can  be  evaluated  approximately  by  discretizing  the  variable 
f and  using  a fast  Fourier  transform  algorithm.  The  spectrum  of  Figure  1(f) 
is  the  result  of  a 256  point  FFT.  We  note  that  both  of  these  estimates 
are  unsatisfactory.*  We  will  see  that  the  new  algorithms  introduced  here 
improve  the  estimated  signal  spectrum  considerably. 


♦Note  that  for  random  signals,  the  periodogram  is  an  inconsistent  estimate. 
Windowing  techniques  may  be  used  to  improve  the  spectrum  estimate  in  the 
sense  that  it  would  be  a consistent  estimate  of  a smoothed  version  of  the 
original  spectrum.  In  this  example,  windowing  did  not  improve  the  situation 
in  so  far  as  the  signal  was  concerned. 


A 
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g.  la:  Spectra  of  Signal,  Clutter  and  Noise 


W 


Fig.  lb:  Actual  Clut1 


17  Samples  of  Clutter  + Signal  + Noise  Fig.  1e:  Max  Entropy  Estimate  Corresponding  to  Id 

(8th  order  model ) 


85.2 


rREO 
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III.  EXTRAPOLATION  OF  BANDLIMITED  SIGNALS 
3.1  Continuous  Time  Signals 

Suppose  we  have  a continuous,  bandlimited  function  f(t)  so  that  its 
Fourier  transform  satisfies 

F(a))  = 0 for  |w|  > a . (15) 

Let  gQ(t)  be  a time  limited  segment  of  f(t)  which  is  available  as  a 
noise-free  observation,  viz., 

'f(t)  , |t|  < T 

go(t)  = • (16) 

0 ■,  |t|  > T 

The  problem  is  to  extrapolate  gQ(t)  outside  the  interval  [-T,T].  This 
is  the  classical  problem  of  extrapolation  of  analytic  functions.  The 
existence  of  a unique  solution  can  be  established  by  observing  that 
bandl imitedness  of  f(t)  implies  it  is  analytic.  This  means  all  its 

derivatives  exist  and  are  bounded  so  that  from  the  Taylor  series  expansion 

.2 

f(T+A)  = f(T)  + Af(T)  + ^ f"(T)  + ...  (17) 

one  can  evaluate  f(t)  outside  [-T,T].  In  practice,  (17)  is  not  very 
useful,  because,  not  only  does  the  series  have  to  be  truncated,  but  also, 
the  evaluation  of  various  derivatives  is  a noise  sensitive  process.  An 
alternative  algorithm  suoqested  by  Slepian,  et  al . [14],  uses  a series 
expansion 

f(t)  = I a„^-(t,Ta)  (18) 

n=0  " " 
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where  {<>^(t,To))  is  a special  set  of  complete  orthonormal  bandlimited 
functions,  called  the  prolate-spheroidal  wave  functions  (PSWF),  which  are 
defined  for  all  t.  The  coefficients  {a^l  can  be  evaluated  as  projections 
of  the  known  function  gQ(t)  on  the  basis  functions  Once  {a^}  are 

known,  the  right  side  of  (18),  considered  valid  for  all  t, gives  the 
extrapolated  signal.  In  practice,  this  method  also  suffers  from  noise 
limitations  and  errors  due  to  truncation  of  the  series.  Moreover  it 
is  extremely  difficult  to  accurately  generate  the  basis  PSW  functions 
so  that  extrapolation  in  a practical  situation  is  quite  hopeless.  For  a 
sinole  example,  see  Frieden  [19]. 

Recently,  Papoulis  [3]  has  introduced  an  iterative  scheme  that  appears 
to  do  better  than  the  PSWF  expansion  method.  The  algorithm  has  the 
following  steps.  The  first  step  is  to  compute  the  Fourier  transform  of 
g^lt)  as  Gg(u))  and  define 


F^  (w) 


Gq(u)  , |a)|  < 0 
0 , jo)!  > o 


Compute  fi(t),  the  Fourier  inverse  of  F^(u)  and  let 


f{t)  . lt|  <T 
f,(t),  |t|  > T 


(19) 


(20) 


Then  compute  G^ (w)  = F[g^(t)]. 

This  is  the  first  step  of  the  iteration, 
function 


At  the  n^^ 


step  form  the 


< o 
> a 


(21) 


* 
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Find  F’\F^(a))]  and  form 


f(t)  . Itl  < T 

fn(t)  . |tl  > T . 


(22) 


Papoulis  has  theoretically  shown  that  f^(t)  converges  to  f(t)  as 
n . If  we  define  a band-limiting  operator  as 

Bf(t)  = f(t)  © ((sinat)/iTt)  (23) 


where  (sinat)/irt)  represents  the  impulse  response  of  a low  pass  filter,  and 
we  define  a time-limiting  operator  as 


Df(t)  = 


f(t)  , \t\  <T 
0 , |tl  > T 


(24) 


D = I - D 


where  I = identity  operator 
then  the  foregoing  algorithm  can  be  written  as  [10] 

g„(t)  . g„(t)  • Hg„.,(t) 

f„(t)  ■ f,(t)  ♦ 

H = DB,  G - BD 


or 


f„(t)  . 
9„(t)  . 


n*l 

k=0 


f,(t)  . G""  = I 


I 

k=0 


go(t) 


(25) 


(26) 


(27) 
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Either  gp(t)  or  ^p(t)  may  be  considered  as  the  extrapolated  signal. 

3.2  Extension  to  Discrete  Time  Signals 

Sabri  and  Steenaart  [10]  have  suggested  a discrete  version  of  this 
algorithm,  as  follows. 

Let  y(k)  be  a discrete,  bandlimited  signal  so  that  its  Fourier  trans- 
form (i.e.,  Z transform  evaluated  on  the  unit  circle)  defined  as 

Y(f)  = i y(k)exp(-j2TTfk)  , (28) 

k = -oo 

satisfies  the  relation* 

Y(f)  = 0 , i > |f|  > o (29) 

We  are  given  a set  of  time  limited,  noise  free  observations 


go(k)  = { 


y(k) 

0 


-M  < k < M 
otherwise  . 


Given  {gjj(k)},  the  problem  is  to  find  an  estimate  of  y(k)  outside  the 
interval  [-M,M].  Following  section  3.1,  we  define  infinite  vectors 

y = [...y(-k)...y(-l),y(0),y(l) y(k)...f 

T (30) 

g„  = [...g„(-k)...g,(-l).g,(0).g„(l).....g,(k)...] 

where  g^  - [0,0.,,,0  gQ(-M) ,gQ(-M+l ) , . . . ,gQ(-l ) .gQ(0) ,gg(l ) t . . .gQ(M) .0 ,0, . . .0] 
We  also  define  a band-limiting  operator  L,  and  a time-limiting  operator  W, 
as  infinite  matrices 


♦This  implies  y(k)  has  been  oversampled  with  respect  to  its  Nyquist  rate. 

This  occurs  quite  often  when  a system  observes  signals  over  a wide  bandwidth. 


A 


I 


(37) 
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In  practice,  the  infinite  matrix  operator  G is  replaced  by  a finite 
matrix,  of  size,  say,  NxN  where  N » (2M+1).  Later  in  this  section  we  will 
show  that  (I-G)  is  singular  for  N = «>,  but  its  finite  approximation  (N<») 
is  non  singular. 


Existence  and  Convergence 

Following  Papoulis,  it  can  be  shown  that  the  above  iterative  algorithm 


satisfies  the  inequality 


/■% 

|Y(f)  - > |Y(f)  - G^(f)|2df  (38) 

-%  -% 


which  says  that  the  mean  square  error  is  reduced  at  each  step.  However,  the 
extrapolated  signal  need  not  converge  to  the  original  signal  y(k)  because 
the  time  limited  discrete  sional  does  not  have  the  analyticity  property  that 
the  continuous  signal  has.  Indeed,  as  we  show  in  the  next  section,  the  above 
(discrete  extrapolation)  algorithm  converges  to  a least  squares,  minimum 
norm  solution  associated  with  the  solution  of  the  equation 


WLy  = gg 

In  terms  of  computational  complexity,  the  iterative  algorithm  requires 
about  4nNlog2N  real  operations  (one  operation  = one  multiplication  and  one 
addition),  where  N is  the  size  of  the  extrapolated  vector  (and  is  much  larger 
than  M)  and  n is  the  number  of  iterations. 

If  the  extrapolation  matrix  E„is  used,  then  once  it  has  been  computed, 
it  requires  j(2M+l ) (N-2M-1 ) operations  to  evaluate  the  extrapolated  signal. 
However  a large  (NxN)  matrix  which  is  ill  conditioned, has  to  be  inverted. 

(see  next  section) 


-22- 
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IV.  EXTRAPOLATION  OF  DISCRETE  TIME,  BANDLIMITED  SIGNALS 

Before  proceeding  to  prove  several  results  related  to  extrapolation 
of  discrete  signals  we  first  consider  several  definitions.  Let  A be  an 
arbitrary  mxn  matrix  and  consider  the  system  of  equations 

Ay  = z (39) 

where  y and  z are  nxl  and  mxl  vectors  respectively. 

4.1  Definitions 

Definition  1:  Least  Squares  Solution 

A least  squares  solution  of  (39),  denoted  by  y,  is  such  that 

llz  - Ay||^  = (z-Ay)^(z-Ay)  (40) 

is  minimum.  This  solution  must  therefore  satisfy  the  equation 

a''’ Ay  = a''’z  (41) 

If  A^A  is  nonsingular  (i.e.  n < m and  rank  of  A is  n)  then 

T T 

y = (a'a)  a z (42) 

is  the  least  squares  solution  and  is  unique.  If  m = n and  A is  nonsingular, 
then 

y = y = A‘''z  (43) 
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If  m < n,  then  A A is  necessarily  singular  and  has  rank  at  most  m.  Then 
(41)  does  not  have  a unique  solution. 

Definition  2:  Minimum  Norm  Least  Squares  Solution 
Let  y^  denote  this  solution.  Then  y^  must  be  that  solution  of  (41) 
which  has  the  minimum  norm  ||y||  . Thus 


y^  = min{||y||^;  A^Ay  = A^z) 


Clearly,  if  rank  A^A  is  n then  y^  = y.  The  minimum  norm  is  simply 
a constraint  that  makes  the  least  squares  solution  unique  for  an  arbitrary  A. 
Definition  3:  Pseudo  Inverse 

We  call  A^  the  pseudo  inverse  of  A[18],  if  for  every  equation  (39), 
the  associated  minimum  norm  least  squares  solution  is  given  by 


y = A z 


This  pseudo  inverse,  also  called  the  generalized  inverse,  satisfies 
the  relations 


When  rank  A = n , 


AA"^  = (AA'*')^ 
A'^A  = (A'^’A)^ 
AA'^A  = A 
A'^AA'^  = A"^ 
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A 


+ 


T T 
(A^A)  a' 


(47) 


If  rank  A = m,  then 


A^(AA^) 


-1 


(48) 


Definition  4:  Singular  Value  Expansion  [21] 

In  general  an  explicit  expression  for  A^,  of  the  type  of  (47)  or  (48) 
is  not  available.  However,  A^  can  be  expressed  as  an  expansion.  Consider 
the  eigenvo^ue  problems 


• V»k 


(49) 


where  k = l,2,...,p  and  p is  the  rank  of  A.  The  vectors  4>|^  and  are  of 
sizes'  nxl  and  mxl  respectively.  Since  A^A  and  AA^  are  non-negative  matrices, 
these  eigen-vectors  exist  and  can  be  orthonormal i zed  so  that 


'^k'^t  ^ '^k,e 


(50) 


From  this,  one  can  express  the  rectangular  matrix  A by  the  expansion,  called 
the  singular  value  expansion,  as 
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= J/kVk 


(51) 


where  > 0,  are  called  the  singular  values  of  A. 
The  pseudo  inverse  can  now  be  written  as 


«+  V 1 T 


(52) 


4 . 2 Properties  of  L 

Now  we  consider  some  useful  properties  of  the  low  pass  filter  operate*- 
defined  in  (31). 

Property  I:  L is  a symmetric  operator,  i.e.,  L = L^.  This  follows 

from  its  definition. 

Property  II:  The  Fourier  spectrum  of  L is  given  by 


e(f)  = 


1 , 0 < |f|  < a 

0,  otherwise. 


(53) 


where  ' J 5 ^ $ 7-  l^^is  is  obvious  since  L is  the  Toeplitz  matrix 
formed  by  the  Fourier  inverse  of  il(f),  the  lowpass  filter  transfer  function. 

Property  III:  Let  S be  a (2M+l)x<»  matrix  operator  whose  elements 
are 


1,  i=j=0±l,t2,...,tM 

0,  otherwise 


(54) 


Basically  S selects  (2M+1)  elements  from  an  infinite  vector.  Consider 
the  (2M+1 )x(2M+l ) matrix 

' T 
L = SLS' 


(55) 


This  is  obvious  because  ideal  low  pass  filtering  a signal  once  is  the 
same  as  doing  it  twice,  i.e., 

Ly  = L{Ly) 

Note  this  implies  the  spectrum  (or  eigenvalues)  of  L must  be  composed 
of  zeros  and  ones  only  [see  (53)]. 


Property  V:  For  every  M < *,  L is  positive  definite.  This  follows 

by  noting 


fO  M 

= I I : 

> m=-M 
-a 


X df 

I ni 


^r|x„(f)pdf.  X„(f)  ^ I 

tn=-M 


> 0,  if  M < 


If  M = “,  then  x^^^  = gives  Xn(|(f)  = 6(f-C)  so  that 


T 

X Lx*  = X Lx*  = 


1,  Ul  < a 

0,  U|  > a 


and  is  not  positive  definite.  Thus,  all  the  eigenvalues  of  L are  positive 


A(L)  >0,  M < oo 


Property  VI:  Let  A (L)  denote  the  largest  eigenvalue  of  L.  Then 

n»x 


A (L)  c 1 , M < 
max'  ' 


1.  M = 


Thus,  for  any  finite  M,  the  eigenvalues  of  L are  bounded  in  the  interval 


(0,1)  i.e.. 


To  prove  (62)  we  note 


0 < X(L)  < 1 , M < 00 


lx)  x'x* 
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From  (59)  we  can  write 


X (L)  = max 
max 


r.j*" 


(f) 


-df 


-df 


(64) 


Since,  for  M < «>,  X|^(f)  is  the  Fourier  spectrum  of  a time  limited  signal, 
it  cannot  be  zero  on  any  finite  interval.  Hence 


•o 

-a 


x„(f) 


V M < “>,  0 


< 


1 

2 


or 


y y X t X < y 

£ i m m-n  n „£■ 


m n 


m=-M 


V M < “>. 


When  M = <»,  one  could  maximize  (64)  by  choosing  a bandlimited  signal  so  that 
the  above  inequality  will  become  an  equality.  This  proves  (62). 

4.3  Iterative  Extrapolation 

With  the  above  definitions  and  properties  we  are  now  ready  to  prove 
the  following  results.  Let  y(k),  k=0,il,...  be  a discrete  time  bandlimited 
signal  as  defined  ir  (28)  and  (29). 

Furthermore,  let  this  signal  be  observed  without  any  noise  over  a 
finite  interval  and  define 


z(k)  = y(k),  - M < k < M 


(65) 


If  z denotes  a (2M+l)xl  vector  and  y is  the  infinite  vector  of  {y(k)},  then 


z « Sy 
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I 


Since  y is  bandlimited,  it  must  satisfy 


Ly  = y 


so  that  we  can  write 


z = SLy 


(66) 


Theorem  1 : Minimum  Norm  Least  Squares  Extrapolation  Theorem 
The  iterative  solution  [see  (32)  to  (37)  and  (57)] 


■q+1  - ^ 

G = L(I-s''’s)  = L(I-W) 


1.2.... 


(67) 


f,  = Lgn  = LS'z 


converges  to  the  minimum  norm  least  squares  solution  y of  (66).  Moreover, 
(67)  is  a special  case  of  a gradient  algorithm  associated  with  the  minimum 
nonn  least  squares  optimization  problem. 

Proof:  An  iterative  gradient  algorithm  associated  with  the  minimum  norm 
least  squares  solution  of  the  general  equation  (39)  is 


>'q.l  ■ ^ J 


. 1 . 1 


(68) 

(69) 


It  is  known  that  y^  converges  to  y as  q under  the  following  conditions  [12]: 
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a) 


(70) 


b)  The  initial  guess  Yq  must  lie  in  the  range  space  of  A^A  e.g.,yQ  = 0. 
From  (66),  letting 


A = SL  (71) 

we  get 

A^A  = lVsL  = LWL  (72) 

A^z  = lVz  = LS^z  = (73) 


Hence  (69)  becomes 


yq+1 


- f,  + y„  - ^ LWL  y„ 
0 1 q 0 g 


(74) 


Now  letting 


yo 


= 0 


and  noting  that  is  bandlimited  i.e.. 


Lf 


1 


(75) 


it  is  easily  verified  by  induction  that  {y^}  is  a bandlimited  sequence  i.e.. 
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Ly^  = y^,  q»1.2.3... 


Using  this  in  (74)  we  get 


^g+l  ” a ^1  ^ ' a '"^^q 


= - f.  + (I  - - LW)y„ 
o 1 o q 


For  0 = 1 , (77)  becomes  the  same  as  (67).  Now  it  remains  to  show  that  this 
algorithm  converges  for  a = 1.  From  (70)  and  (72)  this  requires  us  to  find 
the  largest  eigenvalue  of  LWL.  Now 


‘max"-"-’  ■ • >,»x<'«'> 

■ »™x<'> 


< 1 , V M < ». 


where  we  have  used  Property  VI.  Therefore,  convergence  of  (77)  is  achieved 
whenever 

0 < i < 2 < ^ (78) 

max 

Hence  for  a = 1,  (77)  converges.  This  completes  the  proof  of  Theorem  1. 

An  interesting  question  raised  by  the  foregoing  result  is  "What  is  the 
optimal  value  of  a?"  In  other  words,  we  want  to  find  the  "steepest  descent" 
for  the  gradient  algorithm.  Defining  the  error  vector  at  iteration  step  q as 
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! 

« 


♦ 

I 


e,  ■ y - y,  (79) 

1 

and  noting  that  f-j  can  also  be  written  as 

= Lgp  = LWy  = LWLy  (80) 

we  obtain  from  (74) 

Vi  ■ 

1 q+1 

= (I  - ilWL)  e^.  e^  - y (82) 

This  shows  the  convergence  rate  of  this  extrapolation  algorithm  is  linear. 

Slow  convergence  of  this  algorithm  has  also  been  noted  experimentally 
by  us  and  by  Papoulis  [8]  and  Sabri  et  al  [10].  Since  this  is  a gradient 
algorithm,  convergence  can  be  improved  by  adjusting  a at  every  iteration. 

The  optimal  value  is  given  by 

o„''  = } A = SL.  (83) 

^ h„'A'Ah„ 

q q 

where  h^  is  the  gradient  at  step  q,  defined  as 

hq  = A^(z-Ayq) 

= fl  - LWyq 

This  requires  additional  computations  at  every  iteration  step, 
value  of  a is  desired,  it  is  given  by  [24] 


(84) 

If  a constant 


4 
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= 2/[^n,ax('-WL)  - X^in^LWL)] 


Since 


‘mmC-w-'  ■ » 

‘m.x"-"-!  ' ' 


we  can  take 


r 2 < 2/X 


max 


(85) 


From  our  foregoing  analysis  we  conclude  the  following  about  Papoulis' 
iterative  algorithm. 

1.  The  solution  converges  to  a minimum  norm  least  squares  solution.  Note 
that  continuous  version  of  the  algorithm  converges  to  the  original  band 
limited  signal  y(t),  as  proven  by  Papoulis  [8].  This  reinforces  the  fact 
that  time  limited  discrete  observations  of  a band  limited  signal  need  not 
give  its  exact  extrapolation. 

2.  The  algorithm  is  a gradient  algorithm.  Hence  its  convergence  is 
linear  and  slow.  It  could  be  improved  by  going  to  the  steepest  descent 
algori thm. 

4.4  The  Extrapolation  Matrix 

Now  we  consider  the  extrapolation  matrix  suggested  by  Sabri  and  Steenaart 
[10].  This  is  the  doubly  infinite  matrix  defined  as  [see  eqns.  (32)  to  (36)] 


E = 
00 


(I-G)"'' 


(86) 


r' 


U 

I 

I 

i 

t 


i 

i 

! 

i 


t 


I 

i 
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G = L(I-W) 

In  a practical  situation  the  matrix  6 is  truncated  to  a finite,  but  large, 
NxN  matrix,  say  defined  as 


and  the  corresponding  extrapolation  matrix  is 


We  intend  to  show  that  for  every  finite  N,  exists.  However,  does 
not  exist.  Thus  as  N goes  to  infinity  the  sequence  {Ejij}  becomes  an  ill- 
conditioned  set  of  matrices. 

Lemna  1 : For  every  finite  N,  the  matrix  defined  as 


- I - 


“-n^n 


(88) 


is  nonsingular.  At  N = <»,  is  singular. 

Proof:  From  Property  V,  the  finite  NxN  matrix  is  positive  definite. 
Now  consider  the  symmetric  matrix 


PnS 


Ln  - ^N 


4*  1 W I 

‘-n'^n'-n- 


(89) 


1 

i 

i 

j 


I 


i 


1 

I 


! 


3 


i 


t 

r 


Since  all  the  eigenvalues  of  lie  in  the  interval  (0,1)  we  have 


r 
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. I 

(90) 

1 

Therefore,  for  any  Nxl  vector  x. 

1 

xV^x*  > x\„^X* 

I 

1 

Also 

1 

T M M 

x'l  JmL„x*  = y y X * > 0 

iiNN  m^un-M  nini,nn 

ni=-n  n--n 

where  are  the  (2M+1  )x{2M+l ) elements  of  which  is  positive 

definite.  Clearly,  then  C|^  is  positive  definite.  Hence  exists 

and  is  nonsingular. 

At  N = ®,  L|^is  singular.  Consider  the  eigenvectors  of  the  equation 

LWLx  = Xx 

(91a) 

\ 

« 

Since  LWL  is  symmetric  and  is  of  rank  (2M+1),  there  exists 

an  X such  that 

LWLx  =0,  X ^ 0 

(91b) 

Moreover,  for  every  such  vector  there  exists  a band-limited 

X 

Lx  = X 

1 

< 

which  is  also  a solution  of  (91a).  Thus,  for  all  such  x we 

have 

\ 

i 

4 

* 

-■  ^ 

x^Px  = x^PLx  = x^Lx  - x^L^x  + x^LWLx 
= x^Lx  - x^Lx  + x^LWLx 
= x^LWLx 


Thus  P is  singular. 

4.5  The  Generalized  Inverse 

Having  noted  that  the  foregoing  approaches  give  a minimum  norm  least 
squares  solution,  one  may  attempt  to  find  it  directly.  We  recall  that  the 
given  system  of  equations  is 


SLy  = z 


Defining  A = SL  to  give 


AA^  = SLlV 


= SLS 


We  note  that  L is  positive  definite  (Property  V).  Hence  from  Definition  3 and 
Eqn. (48)  we  can  write  directly  A*  = A^(AA^r^  which  gives  the  extrapolation 


matrix 


E^.  = lV(sls^)’^ 
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, 

= lV(SLS^)'^z  (95) 

This  form  of  solution  was  obtained  for  the  extrapolation  problem  by  Cadzow 
[11]  by  a different  route.  This  method,  we  believe,  is  easier  and  more 
direct.  We  note  that  while  the  extrapolation  matrix  did  not  exist,  the 
extrapolation  matrix  exists.  Moreover,  L is  only  a (2M+1 )x(2M+l ) 

Toeplitz  matrix  so  that  its  dimensionality  is  much  smaller.  However,  as  M 

becomes  large  or  for  certain  combinations  of  M and  o,  L could  be  ill  I 

conditioned.  Experimentally,  the  ill  conditioning  can  be  reduced  by  adding  ! 

a small  diagonal  term  to  L.  This,  however,  will  degrade  the  extrapolated 
estimate. 

4.6  Discrete  Prolate  Spheroidal  Wave  Functions  and  Singular  Value  Expansion 
In  an  earlier  section  we  had  mentioned  that  a continuous  band-limited 
signal  could  be  extrapolated  outside  its  observation  interval,  exactly,  via 
the  PSWF  expansion.  For  the  case  of  discrete  signals, a similar  expansion 
is  possit’e  for  the  minimum  norm  least  squares  extrapolated  estimate.  This 
is  achieved  via  the  singular  value  expansion  described  in  Definition  4. 

Papoulis  and  Bertram  [20]  have  introduced  the  discrete  PSWF  earlier  for 
realization  of  digital  filters  whose  impulse  response  is  an  all-zero  model. 

However,  they  have  not  shown  any  extrapolation  properties  of  these  PSWFs. 

Here  we  introduce  the  discrete  PSWFs  which  also  extrapolate  a discrete 
band  limited  signal  (known  over  a finite  duration)  to  an  infinite,  minimum 
norm  least  squares  signal . 

Following  Definition  4 for  A = SL,  we  consider  the  eigenvalue  problems 


r.>. 


t 


(96) 


= LWU|^  = LS  - M < k < M 


Aa\^  = - M < k < M 


k "k^k’ 


where  > 0,  and  {4>|^}  are  ■»  x 1 and  are  (2M*l)xl  orthogonal  vectors  i .e 


\\  = ^k.e 


From  (96),  mu^  be  a band-limited  signal  satisfying  the  condition 


L(Jk  = \ 


because  = (1/'|^)L(WL  tj^) , and  Lz  is  band  limited, Vz.  Now  define 


Then  (96)  gives 


Ck  = 


^k  = 


SLS'Ck  = \^^k  = \^k 


This  shows  satisfies  the  same  equation  as  Hence  and  must 
be  proportional  i .e. 
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Using  the  orthogonality  condition  (98)  and  the  eqns.  (99),  (100)  and  (97) 
we  find 


or 


(102) 


giving 


Also  from  (96),  this  yields 


\|^  > 0.  - M < k < M 


(103) 


M k < M 


(104) 


Equation  (103)  states  that  the  (2M+l)xl  vector  is  simply  obtained  by  selecting 

the  (2M+1)  elements  {t|^{m),  - M < m < M}  of  and  scaling  then  by 

Eqn.  (104)  is  remarkable  in  that  the  “■xl  vector  is  obtained  by  simply  low 

-1/2 

pass  filtering  the  sequence  {\J»|^(m)}  and  scaling  the  result  by  . This 

means  is  the  extrapolation  of  obtained  by  simple  low  pass  filtering  and 
scaling.  Also  noteworthy  is  the  fact  that  the  sequence  {C|^(n),  - * < m < <*>) 
is  orthogonal  over  the  interval  - M < m < M as  well  as  over  the  infinite  interval. 
This  property  is  similar  to  that  of  the  continuous  PSWFs.  The  (2M+l)x1  vectors 
are  easily  obtained  by  solving  the  eigenvalue  problem  of  (97)  i.e. 
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The  functions  {ij^|^(n)}  and  {(()|^(n)}  have  been  obtained  by  Papoulis  [20]  and  are 

called  the  discrete  PSWF.  However,  their  usefulness  in  obtaining  minimum  norm  i 

least  squares  extrapolation  has  not  been  noted  earlier.  Our  extrapolation  * 

algorithm  requires  the  following  steps.  First  calculate  (2M+1)  orthonormal i zed  1 

vectors  (each  of  size  (2M+1))  by  solving  (105).  Next  obtain  the  (2M+1)  i 

infinite  size  vectors  4)|^  according  to  (104)  by  low  pass  filtering  and  scaling 
of  Then  (108)  and  (109)  give  the  extrapolated  estimate  from  the  observations 
(y (m) , - M < m < M} . 

Properties  of  Discrete  PSWFs:  We  now  summarize  properties  of  these  functions. 

1.  Let  L be  a low  pass  filter  operator  and  L be  a (2M+1 )x(2M+l ) principal 
minor  of  it.  i .e. 

*i.0  ■ s.j  ■ ‘t-J  '-j-o-i'-'-Wl 

Then  the  orthonormal  i zed  eigenvectors  of  L,  denoted  by  form  a complete 
orthonormal  set  of  basis  functions  in  (2M+1 ) dimensional  vector  space.  Let 

" \'*'k  - M < k < M (110) 

X,  >0 

2.  The  discrete  PSWF  are  formed  from  as  ““xl  vectors  {(})|^}  and  are  defined  as 

<{i|^  = 1=  -M<k<M  (111) 

3.  The  PSWFs  are  orthonormal  over  the  infinite  interval,  i.e. 


i 


i 
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4.  The  PSWFs  are  orthogonal  over  the  finite  interval  i.e. 


m=-M  ’ 

5.  The  Fourier  transform  <t>|^{f)  of  4,^(m)}  satisfies  the  eigenvalue  integral 
equation 

-0 

The  proof  is  obtained  by  Fourier  transforming  (111)  and  using  (103). 

6.  The  eigenvalues  of  C lie  in  the  interval  (0,1)  i.e.,  0 < >|^(L)  < 1. 

7.  Let  y(m)  be  a band  limited  signal  whose  spectrum  lies  in  the  interval 
(-o,o).  If  y(k)  is  known  over  the  interval  [-M,M],  then  its  minimum  norm  least 
squares  extrapolation  estimate  is  given  by 

. M a. 

y (iTi)  = I ‘J’k 

k=-M  \ ^ 

(115) 

M 

a.  = y y(m)<i|,(m) 

^ m=-M 

Note  that  this  gives  y^(m)  = y(m)  for  m e [-M,M]. 

So  far  we  have  considered  only  the  case  when  there  is  no  noise  in  the 
observed  signal.  In  the  next  few  sections,  we  consider  other  algorithms  where 
additive  noise  or  interfering  signals  are  allowed  in  the  observations.  The 
algorithms  reduce  to  minimum  norm  least  squares  solutions  in  the  absence  of  any 
noise. 
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V.  A CONJUGATE  GRADIENT  ALGORITHM  FOR  SIGNAL  EXTRAPOLATION 

In  this  section  we  consider  a slightly  different  but  more  general 
problem  than  that  considered  in  the  previous  section.  We  are  interested  in 
extrapolation  and  discrimination  of  two  interfering  signals  (see  Example  1). 

Let  us  assume  that  we  are  given  m observations  which  may  consist  of 
signal,  clutter  and  noise,  where  the  signal  and  clutter  are  bandlimited  in 
mutually  exclusive  bands  (see  Fig.  1(a)).  The  problem  then  is  to  obtain 
estimates  of  extrapolated  signal  and  clutter  outside  their  observation 
interval . 

Let  the  Bandpass  operator  B operate  on  the  signal  in  the  signal  band 
*[f2«^3]  the  low  pass  operator  L operate  on  the  clutter  in  the  clutter 
band  t[0,fi]. 

We  introduce  operator  notation: 
s:  original  signal  (infinite  vector) 

c:  clutter  (infinite  vector) 

n:  noise  (mxl  vector) 

y:  observed  samples  (mxl  vector) 

B:  bandpass  operator  (infinite  matrix) 

L:  Low  pass  operator  (infinite  matrix)  ^ 

S:  Selection  operator  (mx®) 

The  matrices  B and  L are  Toeplitz  matrices  determined  from  the 
sequences  (bj^},  respectively,  as  the  Fourier  inverses  [see  Fig.  11], 

/s  sin(2TTf.,k)  - sin(2nfpk) 

= 2]  cos(2nfk)df  = ^ ^ (116) 

^2  < 
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1.-2  cos  (2Tifk)df 

Jn 


sin{2Tif,k) 


L . 

The  selection  matrix  S,  (introduced  earlier  also)  is  defined  as 

S = [0  ; : 0]  (11! 

where  I is  an  m x m identity  matrix.  The  matrix  selects  the  observed  m 
m 

samples  from  an  infinite  size  vector.  The  observation  vector  can  be 


written  as 


y = Ss  + Sc  + n 


Since  s and  c are  bandlimited  signals  they  satisfy  Bs  - s and  Lc  - c. 


Thus,  we  can  write 


y = SBs  + SLc  + h 


y = S[B  L]  fs'l  + n 


^ Hx  + n 


when  H = S[B  L]  and  x = 

LcJ 

Now  the  problem  is  to  find  an  estimate  of  x,  given  y.  In  the 
absence  of  clutter  (c=0),  the  problem  reduces  to  the  extrapolation  problem 
considered  earlier.  The  solution  k obtained  by  minimizing  the  least 


squares  norm 


J = lly-Hx|l 


is  given  by 
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X = 

where  (H^H)  is  the  pseudo  inverse  operator  of  H provided  exists. 

However,  this  is  not  the  case  as  explained  below. 

The  infinite  matrices  B and  L are  Toeplitz  and  are  diagonalized  by  the 
Fourier  transform  operator.  Therefore, 

F[B  = A{f) 

where  F is  the  Fourier  transform  operator  and  F^  is  its  conjugate  transpose 
and  A(f)  is  as  shown  in  figure  11. 

This  also  implies  that  the  operators  B and  L are  idempotent  i.e., 

2 2 

B = B and  L = L.  From  Fig.  11  the  operator  [B  L]  is  singular  which  in 


turn  implies 


T 

h'h  = ' 


S'S[B  L] 


(123) 


is  singular,  and  has  a rank  of  atmost  m, 

A number  of  approaches  are  once  again  possible  to  find  the  pseudo 
inverse  of  H,  as  discussed  in  the  previous  section.  Here  we  consider  a two 
step  gradient  method. 

An  example  of  a two  step  gradient  method  where  initial  convergence  is 
extremely  rapid  is  the  Conjugate  Gradient  Method  [12,13]  which  is  based 
on  the  following  ideas.  Let  Q denote  the  matrix  H^H  (which  has  dimension 
2Nx2N  to  extrapolate  the  signal  and  clutter  each  to  N points). 

Let  t = 2N.  A set  of  I vectors  {d.}  are  conjugate  if 


d^^Qdj  = 0 , i / j . 


A 


EX. 
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/ 


Since  Q is  symmetric  a set  of  such  vectors  exist  and  forms  a basis.  The 
solution  vector  can,  therefore,  be  written  as 


' ■ ■ 


The  scalers  {a^}  and  vectors  {d^}  must  be  found  in  a computationally 
feasible  way.  One  way  is  via  the  algorithm  [13] 


-k+l  -k  ^ . 

X = X + a^d|^ 


(124) 


“iT’wk 


where  the  {dj^}  are  generated  by 


<^k.l  = -'k.l  " ¥k 


(125) 


®k  = fT 


^k+1  ‘^‘^k 


'k  Q^k 


The  vectors  {Cj^}  are  the  gradients  of  J at  each  iteration 


+ Vl^Vl 


(126) 


and  the  initial  conditions  are 


-1 


X = y and  d^  = C^. 


(127) 


The  minimum  is  achieved  in  at  most  8.  steps,  and  the  method  is  step  for 
step  better  than  a gradient  method.  A very  attractive  feature  of  this 


algorithm  is  that  large  reductions  in  error  are  achieved  in  the 
first  few  steps  [13]. 

Looking  at  the  computations  involved,  we  see  that  except  for  the  single 
matrix  vector  product  Qdj^,  all  vector  operations  involve  only  i multiplica- 
tions. Hence,  the  order  of  computations  for  each  iteration  is  4mN. 

(Note  that  Q is  composed  from  Toeplitz  matrices,  so  that  advantage  of  FFT 
method  could  be  taken  to  evaluate  Qdj^.)  When  Q is  highly  ill-conditioned, 
or  singular,  the  iterations  must  be  stopped  at  an  optimum  point.  Alter- 
natively, we  may  add  a small  value  c of  the  order  of  10  ^ to  each  diagonal 
term  of  Q.  This  minimizes  convergence  ambiguity  due  to  ill-conditioning 
and  stabilizes  the  iterations  greatly. 


VI.  A MEAN  SQUARE  EXTRAPOLATING  FILTER 


With  the  formulation  and  notation  of  the  problem  of  Section  V,  we 
have^ 

y = Hx  + n . (128) 

Now  we  assume  that  x is  a random  Gaussian  vector  whose  autocorrelation 
matrix  is  denoted  by  R^.  The  minimum  mean  square  estimate  of  x is  given 
by  the  Wiener  filter  as* 

X = [E(x  y^)  ][E(y  y^)]'\  ^ Gy  . (129) 


Assuming  noise  to  be  independent  of  x,  it  is  easy  to  obtain 

G = Ry(HRy+R^)‘^ 

= E(nn^ 

which  is  equivalent  to  the  equations 

(HR  H^+R  )z  = y 
X n 

= ECxx"^]  = E[(^)(sV)] 

E(ss^)  E(sc^) 

E(cs^)  E(cc^) 


(130) 


• (140) 


* We  note  here  the  similarity  between  the  extrapolation  problem  and  the 
restoration  problem  in  image  processing  [17],  For  example,  an  image 
blurred  by  a low  pass  type  operator  H and  contaminated  by  additive  noise 
would  give  rise  to  an  equation  similar  to  (128). 

* Here  E is  the  expectation  operator. 
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If  signal  and  clutter  are  uncorrelated  (non  overlapping  power  spectra) 


E(ss')  0 


E(cc^) 


In  the  special  case  when  the  noise  and  clutter  are  absent,  we  have 


H = B 

Rn  = 0. 

Thus,  we  have  z = V-  Since  s(k)  is  a bandlimited 

random  process,  we  must  have 


which  gives 


BR^B  -R^ 


X = R^s''’(SR^s'^)'''y 


In  the  worst  case,  when  we  do  not  know  R^,  we  can  simply  assume  the 
power  spectrum  of  x(k)  to  be  flat  in  its  bandwidth  so  that 


Rx  = B 


and  9.  = B5'(5BS  )^y  where  (SBS')^  is  the  nseudo  inverse  of  (SBS^).  This  is 
the  same  result  as  obtained  by  Cadzow  [See  Eqn.  (95)].  Thus  Cadzow's  one 
shot  method  is  a special  case  of  this  extrapolation  filter.  In  the 
presence  of  noise,  the  extrapolation  filter  estimate  is 

X = BS^(SBS^+R^)‘^y 

where  (SBS^+R  ) ^ exists  and  is  unique, 
n 
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VII.  A RECURSIVE  EXTRAPOLATION  ALGORITHM 

In  this  section  we  present  a recursive  least  squares  algorithm  based 
on  Kalman  filtering  techniques  where  the  extrapolated  signal  estimate 
is  updated  recursively  as  a new  observation  sample  arrives. 

Based  on  the  formulation  of  the  problem  as  in  section  V we  rewrite 
equation  (121); 


y = Hx  + n 

The  r nervation  yj^  can  be  written  as 

^k  ~ ^k^  ^ ^k  k=0,2 , . . . ,m-l  (142) 

where  h^  is  the  k^^  row  of  H and  n.  is  zero  mean  white  Gaussian  noise, 
k 

The  state  equation  for  the  unknown  extrapolated  vector  x can  be  written 
as 


^k+1  " ’^k 


(143) 


with  initial  condition  x^  = x,  where  x is  a random  vector  whose  covariance 
is  given  by 


P = cov(x^)  = cov(x)  = H = (a 
0 0 ' 


k-e,' 


(144) 


Since  H is  idempotent  i.e.  HH  = H,  it  is  easy  to  verify  that 
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(145a) 


(145b) 


(145c) 


The  Kalman  filter  associated  with  equations  (142)-(144)  is  the  recursive 
least  squares  filter  (e.g.  see  Nahi  [16]) 


Vi  = ^ = 0 


where  is  the  k estimate  of  x and  is  the  Kalman  filter  gain. 
The  associated  Riccati  equation  is 


Vl  = - E7  - ^7  Vk^k)  " Vn^C^)^  (147) 


g,.  = — P,  h, 
k c,  k k 
k 


Equation  (147)  then  reduces  to 
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Assuming  that  we  are  basing  our  extrapolation  on  m observed  sample  values  and 

that  we  are  extrapolating  to  N points,  the  maximum  storage  required  is  mN. 

The  major  computation  is  in  the  calculation  of  the  (153)  and  it  is 

2 

easily  seen  that  the  order  of  computations  involved  is  '0(m  N). 
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VIII.  EXAMPLES,  RESULTS  AND  COMPARISONS 

We  will  now  consider  several  examples  to  study  the  performance  of 
the  foregoing  algorithms.  In  Example  1 below,  we  consider  the  observations 
to  be  of  the  form 

y(k)  = s(k)  + c(k)  + n(k)  , -M  ^ k jc  M 

where  c(k)  is  the  bandlimited  clutter  sequence,  and  s(k)  is  a bandlimited 
signal  whose  band  limits  are  known,  i.e. , 

S((i))  = 0,  (1)  ,1^2]. 

2 

n(k)  is  a zero  mean  white  noise  process  with  variance  . 

In  examples  2-9  the  observations  are  of  the  form 

y(k)  = s(k)  + n(k)  ; -Mj<kj<M  . 

The  following  examples  have  been  considered. 

la.  y(k)  = s{k)  + c(k)  + n(k)  ; -8  ^ k 8 

s(k)  = 1 . 69sin( . 3971k) 

= 0.38271 
11)2  = 0.  39771 

clutterband:  = [0,0.167t] 

% ■ “ '3 

SCR  = -4.1dB,  SNR  = 19dB 

lb.  y(k)  = s(k)  + c(k);  - 8 < k < 8 

W-j  = .377 


0)2=  .477 
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0)^  = .875^2 
s(k)  = 1 .69sin(a)jk) 


2.  s(k)  = sin(.99u^k)  + sin(.85oj2k) 


■ u)^  = 0.8co2 

ii>2  - 2ti/50 

3.  s{k)  = sin(.99u)2k)  + sin(.85u)2k) 


cj-i  = 0 

- 0-1^ 

The  above  two  cases  have  the  same  signal  , but  the  knowledge  of  band- 
limits  is  different.  In  the  following  cases  we  have  other  signals  with 
different  bandlimits.  In  examples  8 and  9 we  have  additive  noise  also. 

4.  s(k)  = sin(.365TTk)  + sin(.385TTk) 


= .371 
W2  = -477 


5.  s(k)  = sin(.365iTk)  + sin(.3857ik) 


6.  s(k)  = sin(.985oj2k)  + sin{.975a)2k) 

~ 0.  96u)2 

= 0.4it 

o 2 = 0 
n 

sinu)„k 

u)l  - 0 

^2  ~ •If' 

= 0 

sino)  k 

8- 

= 0 

t02  = 0.  Iti 

E(n(k))  = 0.  Op^  = 0.01 
SNR  = 7.4dB 

9)  s(k)  = sin(.99w2k)  + sin(.85a)2*^) 

uj^  “ . 8co2 

u)^  = 2tt/50 
E(n(k))  = 0 

o ^ = 0.1  , SNR  = 21.6dB 
n 


In  example  la,  application  of  an  eighth  order  autoregressive  model  to 
the  given  data  yields  the  spectrum  shown  in  Fig  1(e).  The  estimate  seems 


J 


to  be  completely  dominated  by  the  clutter.  A 256-point  DFT-spectrum 
estimate  of  the  data  is  shown  in  Fig.  1(f),  and  clearly,  this  also  is 
most  unsatisfactory.  Figures  1(g)  and  1(h)  show  the  signal  extrapolated 
to  199  points  by  the  Conjugate  Gradient  Algorithm  (using  only  10  itera- 
tions) and  the  Mean  Square  Extrapolation  Filter.  Figures  l(i)  and  l(j) 
show  the  corresponding  results  for  the  extrapolated  clutter.  The  Maximum 
Entropy  method  is  applied  to  the  extrapolated  signal  and  extrapolated 
clutter  separately  using  a fifteenth  order  model.  Fig.  l(k)  and  l(il)  show  the 
extrapolated  signal  spectra  using  the  Conjugate  Gradient  Algorithm  and  the 
M.S.  Extrapolation  filter.  They  yield  the  signal  peak  at  the  correct 
location.  Figures  l(m)  and  l(n)  show  the  corresponding  clutter  spectra. 

In  example  lb,  the  signal  bandwidth  is  increased  from  the  previous  case. 

The  actual  spectra  of  signal  and  clutter  are  shown  in  fig.  2(a).  The  FFT 
spectrum  using  256  points  is  shown  in  fig.  2(b),  and  the  Maximum  Entropy 
Spectrum  using  an  eight  order  model  in  fig.  2(c).  After  extrapolating  signal 
and  clutter  to  125  points  each,  using  the  conjugate  gradient  algorithm  their 
respective  spectra  are  calculated  by  the  Maximum  Entropy  Method  using  15th  and  12th 
order  models  as  shown  in  figures  2(d)  and  2(e)  respectively.  Figure  2(d)  shows 
a peak  at  the  signal  frequency  along  with  two  subsidiary  peaks  which  are  seen 


to  occur  exactly  at  the  filter  cut-off  frequencies.  Figure  2(f)  shows  the 
FFT  spectrum  of  the  extrapolated  signal.  We  note  that  this  is  a much  better 


Algorithm 


fNTROP 


Fig.  Ik:  Max  Entropy  Spectrum  of  (Ig)  (15th  order  model)  Fig,  11:  Max  Entropy  Soectrun  of  (Ih)  (15th  order  model) 


Fig.  Im:  Max  Entropy  Spectrum  of  (li)  (15th  order  Model)  Fig.  In:  Max  Entropy  Spectrum  of  (Ij)  (15th  order  model) 


S 
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The  observed  data  consists  of  17  samples  and  it  is  extrapolated  to 
199  points  using 

1.  Papoulis'  iterative  algorithm 

2.  the  matrix 

3.  the  matrix  E after  it  has  been  stabilized 

c 

4.  the  conjugate  gradient  algorithm 

5.  the  M.S.  extrapolation  filter  for  Examples  8 and  9 when  the  signal 
contains  noise. 

Extrapolation  using  Papoulis'  algorithm  was  done  with  30  iterations  and  with 
10  iterations  using  the  conjugate  gradient  algorithm. 

Since  all  the  algorithms  ultimately  yield  a minimum  norm  least-squares 
solution,  we  expect  equivalent  results  from  all  the  algorithms.  Figure  3(c) 
shows  the  extrapolation  using  Papoulis'  algorithm,  and  is  seen  to  give  a 
reasonable  result.  Figure  3(d)  shows  the  extrapolation  via  the  matrix  E^. 
After  a stabilizing  diagonal  term  of  the  order  of  10  ^ is  added  to  the  matrix 
(5L5^),  the  extrapolation  obtained  is  as  shown  in  figure  3(e).  A close 
examination  of  figures  3(a),  3(d)  and  3(e)  shows  a slight  phase  shift 
of  the  signal  3(d).  This  may  be  attributable  to  the  ill-conditioning 
of  (SL5^).  Figure  3(f)  shows  the  extrapolation  by  the  Conjugate  Gradient 
method  after  10  iterations.  This  algorithm  gives  a slightly  inferior 
extrapolation  compared  to  the  others  in  this  case,  and  perhaps  needs  a few 
more  iterations.  In  later  examples,  however,  it  is  seen  that  10  iterations 
are  sufficient. 

In  the  set  of  figures  4,  we  have  the  same  signal  and  observations,  but 
the  uncertainty  in  bandwidth  is  increased  over  the  previous  case.  The 
relatively  degraded  extrapolations  achieved  now  show  that  there  is  a 
trade-off  between  bandwidth  uncertainty  and  extrapolation  length,  which  is 
to  be  expected.  The  set  of  figures  5,  6,  7 lead  to  similar  conclusions  for 
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signals  with  differing  frequencies  and  different  bandwidths. 

In  the  set  of  figures  8,  we  consider  a sin  x/x  type  of  signal.  The 
phase  distortions  of  the  extrapolated  estimate  in  figure  8(d)  dramatically 
illustrate  the  ill-conditioned  nature  of  (SLS^).  In  the  above  examples 
there  was  no  noise  in  the  observations.  Now,  we  consider  (examples  9 and 
10)  signals  in  additive  zero  mean,  white  Gaussian  noise  at  the  SNR  7.4dB, 
21.6dB  respectively.  The  results  show  that  the  Mean  Square  Extrapolation 
filter,  which  takes  noise  statistics  into  account  produces  the  best 
extrapolation  (Figs.  9(g)  and  10(g))  among  the  algorithms  considered 
experimentally.  Similar  results  are  to  be  expected  form  the  recursive 
extrapolation  algorithm  of  section  VII.  Further  simulations  are  required 
to  study  the  numerical  behavior  of  this  algorithm. 


Given  Observations  fl7  Samples) 


Fig.  3e:  Signal  Extrapolated  After  Adding  a Stabilizing  Fig.  3f:  Signal  Extrapolated  by  Conjugate 

Diagonal  Term  to  Matrix  E Gradient  Algorithm 


e ti’l  3 3ICNt,  i <TR^'-Sl  (,’l  [; 


Fig.  4c;  Signal  Extrapolated  by  Papoulis'  Iterative  Fig.  4d:  Signal  Extrapolated  Via  Matrix 

Algorithm 


trapol ated 


Fig.  5e:  Signal  Extrapolated  After  Adding  a ‘'stabilizing  Fig.  5f:  Signal  Extrapolated  by  Conjugate 
Diagonal  Term  to  Matrix  E Gradient  Algorithm 


Signal  Extrapolated  by  Papoulis'  Iterative  Algorithm  Fig.  3d:  Signal  Extrapolated  Via  Matrix 


Signal  Extrapolated  by  Papoulis'  Iterative  Algorithm  Fig.  7d:  Signal  Extrapolated  Via  Matrix 


Fig.  8a:  Original  Signal  Fig.  8b:  Given  Observations  (17  Samples) 


8c:  Signal  Extrapolated  by  Papoul 
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Fig.  9c:  Signal  Extrapolated  by  Papoulis'  Iterative  Algorithm  Fig.  9d:  Signal  Extrapolated  Via  Matrix 
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Fig.  10a:  Original  Signal  Fig.  10b:  Given  Observations  (17  Samples) 
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IX.  CONCLUSIONS 

The  problem  considered  here  was  to  effectively  discriminate  the 
signal  from  the  interfering  clutter,  based  on  a small  number  of  observation 
samples.  Conventional  techniques  like  the  OFT  and  Maximum  Entropy  Methods 
are  seen  to  yield  poor  results. 

Given  that  the  signal  and  clutter  are  band  limited  in  mutually  exclusive 
bands,  a good  method  oF  improving  resolution  is  to  extrapolate  the  signal 
outside  the  observation  interval  and  then  estimate  the  spectrum.  Though 
continuous  time  band  limited  signals  can  be  extrapolated  exactly  outside  the 
observation  interval,  this  is  not  possible  in  the  discrete  time  case.  In 
fact,  in  this  report  we  have  proved  that  the  discrete  extension  of  the 
continuous  algorithms  leads  to  an  extrapolated  signal  which  is  optimum 
in  a minimum  norm  least  squares  sense. 

We  have  introduced  some  new  algorithms  for  signal  extrapolation,  viz.; 

a)  Conjugate  Gradient  Algorithm 

b)  Mean  Square  Extrapolation  Filter 

c)  Recursive  Extrapolation  Filter 

d)  An  Extrapolation  Algorithm  via  Discrete  PSWFs. 

The  Conjugate  Gradient  method  is  an  iterative  technique  which  has 
a rapid  initial  rate  of  convergence  and  hence  has  an  obvious  advantage 
over  Papoulis'  algorithm  which  has  a linear  rate  of  convergence.  The  Mean 
Square  extrapolation  filter  is  a non-iterative  method  that  takes  noise 
statistics  into  account  and  simultaneously  filters  the  clutter  from  the 
signal.  Cadzow's  one  shot  method  is  seen  to  be  a special  case  of  this 


filter. 


Further  experiments  have  to  be  performed  for  the  Recursive  Extrapolation 
filter  (which  also  considers  noise  statistics)  and  for  extrapolation  via 
Discrete  PSWFs.  In  practice,  extrapolation  is  achieved  only  to  a limited 
extent  beyond  the  observation  interval  and  depends  on  the  signal  bandwidth 
uncertainty.  The  larger  this  bandwidth  uncertainty,  the  smaller  the  length 
to  which  the  signal  can  be  extrapolated.  Although,  in  the  absense  of  noise, 
all  algorithms  yield  a minimum  norm  least  squares  solution,  their  implementations 
are  different  and  lead  to  different  truncation  errors.  Such  error  analysis 
and  other  numerical  features  for  the  extrapolation  problem  need  further 
work.  Our  experiments  show  that  whatever  the  uncertainty  of  the  signal 
bandwidth  might  be,  extrapolation  of  the  signal  followed  by  a spectral  estimation 
improves  the  spectral  estimate. 
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Table  1 


1 

METHOD 

COMPUTATIONAL 

COMPLEXITY 

GENERAL 

COMMENTS 

1 

f 

f 

i 

MAX.  ENTROPY  or 
Autoregressive  (Burg, 
Parzen  & others) 

pxp  Toeplitz  Eqns. 

+ FFT  operation 
- 3p2  + O(NlogN); 
p«N 

1.  Simple  linec-*^  Eqns.,  easy  to 
implement. 

2.  Good  results  for  all-pole  spectra. 

3.  Fails  in  the  presence  of  noise  and 
clutter-unless  a large  number  of 
observation  samples  are  available. 

4.  Order  of  the  model  difficult  to 
select. 

5.  Performance  improved  by  applying  it 
on  extrapolated  signal. 

i 

2 

Continuous  PSWF 
(Slepian  et  al .) 

Very  large.  Func- 
tions extremely 
difficult  to  cal- 
culate 

1.  For  extrapolation  of  bandlimited, 
continuous  signals;  existence  guaran- 
teed. 

2.  Extremely  difficult  to  implement. 

3.  Noise  sensitive 

4.  Useful  in  establishing  existence, 
uniqueness  & other  properties. 

3 

Iterative 

Extrapolation 

(Papoulis) 

0(4Nlog2N)  real 
operations  per 
iteration 

1.  Easy  to  implement 

2.  Is  a gradient  algorithm  with 

linear  convergence.  Requires  a large 
number  of  iterations  and  FFT 
operations  at  each  iteration. 

3.  Does  not  take  into  account  noise 
statistics . 

' 4 

Extrapolation 

Matrix,  Ea> 

(Sabri  and 

Steenaart) 

If  observed  data  = 

2 M+1  and  extrapo- 
lated length  = N 
then  ■■  0(n3)  to  in- 
vert NxN  Matrix 
+1/2(2M+1)(N-2M-1) 
operations  to  find 
extrapolated  signal 

1.  Eoo  does  not  exist.  Em  exists. 

2.  A large  (NxN)  ill-conditioned 

matrix  has  to  be  inverted.  1 

3.  Noise  sensitive.  Can  be  stabilized  by  J 

adding  a diagonal  term  to  G.  j 

4.  Noise  statistics  not  considered.  j 

I 

5 

E^  (Cadzow) 

If  observed  data  = 
m pts  extrapolated 
to  N pts, 
0(3m^+m^+mN) 

1.  Easy  to  implement,  if  Ej-  is  stable. 

2.  An  ill-conditioned  matrix  (mxm  i 

Toeplitz)  has  to  be  inverted. 

3.  Noise  sensitive.  Can  be  stabilized. 

4.  Noise  statistics  not  considered. 

1 

6 

Conjugate 

Gradient  (Jain  & 
Ranganath) 

- 0(2inN)  operations 
per  iteration 

1.  Easy  to  implement 

2.  Extremely  rapid  initial  rate  of 
convergence 

3.  Small  number  of  iteration  required 
in  practice 

4.  Noise  sensitive,  but  can  be  stabilized 

5.  Noise  statistics  not  considered 

' ■ - 1 
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Table  1,  continued 


METHOD 

COMPUTATIONAL 

COMPLEXITY 

GENERAL 

COMMENTS 

7 

Mean  Square 
Extrapolation 

Filter  (Jain  & 1 

Ranganath)  1 

~ t3ni^+m^+mN ) 
operations.  mN 
operation  once 
Filter  gain  has 
been  computed. 

1.  Easy  to  implement 

2.  An  mxm  Toeplitz  matrix  has  to  be 
inverted. 

3.  Takes  Noise  statistics  into  account 

4.  Reduces  to  Ej.  in  the  noise  free  case. 

8 

1 

Recursive  Extra- 
polation (Jain  & 
Ranganath) 

N operations 
per  data  point, 
if  gains  are  pre- 
computed. 

1.  Easily  implemented,  updates  extra- 
polated estimate  as  new  data  arrives. 

2.  Takes  noise  statistics  into  account. 

9 

Discrete  PSWF 

Singular  Value 
Expansion 
(Papoulis,  Jain  & 
Ranganath) 

Requires  solving 
for  eigenvalues 
and  eigenvectors 
of  a (2M+l)x(2M+l) 
Toeplitz  matrix,  a 

1 low  pass  filtering 

1 operation  and  a 

1 finite  series 

1 expansion. 

1.  More  difficult  to  implement  than 
(7)  or  (8).  Easy  to  implement  once 
the  eigenvectors  have  been  computed 

2.  Noise  statistics  not  considered. 

3.  Accuracy  depends  on  the  accuracy  of 
eigenvalues  and  eigenvectors. 

I 

I 

i 


i 


4 
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